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ABSTRACT 

We investigate the time evolution of luminous accretion disks around black holes, conducting the 
two-dimensional radiation-hydrodynamic simulations. We adopt the a prescription for the viscosity. 

The radial-azimuthal component of viscous stress tensor is assumed to be proportional to the total 
pressure in the optically thick region, while the gas pressure in the optically thin regime. The viscosity 
parameter, a, is taken to be 0.1. We find the limit-cycle variation in luminosity between high and low 
states. When we set the mass input rate from the outer disk boundary to be 100 Le/c^, the luminosity 
suddenly rises from O.SLe to 2 Le, where Le is the Eddington luminosity. It decays after retaining high 
value for about 40 s. Our numerical results can explain the variation amplitude and duration of the 
recurrent outbursts observed in microquasar, GRS 1915-1-105. We show that the multi-dimensional 
effects play an important role in the high-luminosity state. In this state, the outflow is driven by the 
strong radiation force, and some part of radiation energy dissipated inside the disk is swallowed by 
the black hole due to the photon-trapping effects. This trapped luminosity is comparable to the disk 
luminosity. We also calculate two more cases: one with a much larger accretion rate than the critical 
value for the instability and the other with the viscous stress tensor being proportional to the gas 
pressure only even when the radiation pressure is dominant. We find no quasi-periodic light variations 
in these cases. This confirms that the limit-cycle behavior found in the simulations is caused by the 
disk instability. 

Subject headings: accretion: accretion disks — black hole physics — hydrodynamics — radiative 
transfer — stars: individual (GRS 1915-1-105) 


1. INTRODUCTION 

Microquasars in our Galaxy display large flux vari¬ 
ability in X-ray band. Such luminosity variations are 
thought to reflect violent phenomena in accretion disks 
around black holes. The dwarf-nova type disk-instability 
is responsible for the luminosity changes on the time- 
scale of months (e.g., Mineshige & Wheeler 1989; for a 
review Kato et al. 1998, §5). In contrast, the mecha¬ 
nism causing the short-term variability (10 — 100 s) is 
not well understood yet. One of the plausible mecha¬ 
nisms for it is thermal and secular instability which arises 
when the radiation pressure becomes dominant over gas 
pressure (e.g., Lightman & Eardley 1974; Shibazaki & 
Hbshi 1975; Pringle 1976; Shakura & Sunyaev 1976). S- 
shaped sequence on the Mace ~ U (mass accretion rate 
vs. surface density) plane was completed by the finding 
of a slim-disk branch (Abramowicz et al. 1998), in which 
the viscous heating is balanced by the advective cooling. 
It was suggested that the disk may undergo limit-cycle 
oscillations, like the dwarf nova outbursts (for a review 
Kato et al. 1998, §10). This instability is expected to 
occur when the mass accretion rate is comparable to or 
moderately exceeds the critical value, Le/c^, where Le 
is the Eddington luminosity and c is the velocity of light. 

The limit-cycle oscillations have been investigated 
by one-dimensional (ID), vertically-integrated approach. 
The quasi-periodic outbursts were first demonstrated by 
Honma et al. (1991) using the time-dependent ID simu¬ 
lations of the accretion disks and have then been investi¬ 
gated in detail (Szuszkiewicz & Miller 1997, 1998, 2001; 
Watarai & Mineshige 2003). However, the ID model can¬ 
not treat the multi-dimensional motion, i.e., convection. 


circulation, and outflow, although they would influence 
the flow dynamics significantly via the transport of mass, 
momentum/angular momentum, and energy. The modi¬ 
fication of the accretion disk model, where the mass ejec¬ 
tion from the disk surface is taken into account, was sug¬ 
gested by Nayakshin et al. (2000) (see also Janiuk et al. 
2000; 2002 and Janiuk & Czerny 2005). But they still 
used simplihed ID method without treating the multi¬ 
dimensional flow motion. We need to perform at least 
two-dimensional (2D) analysis in order to understand the 
accretion flow correctly. 

Recently, Teresi et al. (2004a, 2004b) first repro¬ 
duced the quasi-periodic luminosity variations by a 
2D smoothed particle hydrodynamic (SPH) simulations. 
However, their simulations focused on the time evolu¬ 
tion of the disk itself and do not treat the optically thin 
regime correctly, i.e., the disk surface as well as the at¬ 
mosphere surrounding the disk. They assumed the equi¬ 
librium between gas and radiation without treating the 
energy of gas and radiation separately. In their simu¬ 
lations, the transport of the radiation energy was not 
solved in the optically thin regime, though the radiation 
flux was evaluated by using the diffusion approximation 
deep inside the disk. The thermal energy was extracted 
from the SPH particles located near the disk surface. It 
should be calculated by solving the interaction between 
gas and radiation as well as radiative transfer. Their 
assumptions and method are valid only in the optically 
thick regime, i.e., deep inside the disk. As a result, the 
outflow was not found in their simulations even though 
the luminosity exceeded the Eddington luminosity in the 
high-luminosity state. The outflow would have an im- 
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portant role on the evolution of the disk through the 
extraction of the mass, momentum/angular momentum, 
and energy from the disk. In addition, the simulations 
for stable disks were not performed in their work. A com¬ 
parison of their results with those of simulations for sta¬ 
ble disks is essential in order to understand the physical 
mechanism of the disk oscillations. Hence, it is impor¬ 
tant to confirm the limit-cycle behavior by the grid-based 
simulations, in which the multi-dimensional radiation hy¬ 
drodynamic (RHD) equations are properly solved. It is 
also important to investigate the time evolution of not 
only an unstable disk but also a stable disk. 

The photon trapping is also basically a multi¬ 
dimensional effect, by which some or large part of pho¬ 
tons generated inside the disk is swallowed by the black 
hole in supercritical accretion regime, leading to a re¬ 
duction of the energy conversion efficiency (Ohsuga et 
al. 2002, 2003, 2005). It has been reported by the anal¬ 
ysis of the X-ray data of GRS I9I5-I-I05 that the disk 
luminosity is comparable to or exceeds the Eddington 
luminosity in the high-luminosity state (Yamaoka et al. 
2001), implying the supercritical accretion. Thus, the 
photon trapping is expected to appear in high-luminosity 
state of the quasi-periodic oscillations. 

Here, by solving full set of 2D RHD equations includ¬ 
ing viscosity term, we report the 2D RHD model for the 
quasi-periodic oscillations of the accretion disks around 
black holes. Our numerical simulations carefully treat 
the 2D effects, including the outflow motion and the pho¬ 
ton trapping. We also show in our simulations that the 
disk is unstable on condition that the mass accretion rate 
is moderately larger than the critical value and the vis¬ 
cous stress tensor is proportional to the total pressure. 
It is consistent with the disk theory and proves that the 
bursting behavior in our simulations arises from the disk 
instability in the radiation-pressure dominant region. In 
§2, our model and numerical method are described. We 
present the numerical results in §3. §4 and §5 are devoted 
to discussion and conclusions. 


2. MODEL AND NUMERICAL METHOD 

Basic equations and our numerical method are de¬ 
scribed in detail in Ohsuga et al. (2005). Here we briefly 
summarize the model and the numerical method. The 
set of RHD equations including the viscosity term are 
solved by an explicit-implicit finite difference scheme on 
the Eulerian grids. We use spherical polar coordinates 
(r, 6 , if), where r is the radial distance, 0 is the polar an¬ 
gle, and f is the azimuthal angle. In the present study, 
we assume the axisymmetry as well as the reflection sym¬ 
metry relative to the equatorial plane. We also assume 
that only r^s-component of the viscous stress tensor, r^,^, 
is non zero. This component plays an important role for 
the transport of the angular momentum and heating of 
the disk plasma. We describe the gravitational field of 
the black hole in terms of pseudo-Newtonian hydrody¬ 
namics (Paczynski & Wiita 1980), and assume the flow 
to be non self-gravitating. The basic equations are the 
continuity equation. 
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and the energy equation of the radiation. 
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Here, p is the gas mass density, v = {vr,V 0 ,Vp) is the 
velocity, Pg^s is the gas pressure, e is the internal energy 
density of the gas, B is the blackbody intensity, Eq is the 
radiation energy density, Fq = (Fq , Fg ) is the radiation 
flux, Po is the radiation pressure tensor, rj is the dy¬ 
namical viscosity coefficient, k is the absorption opacity, 
and x{= K + paTlrrip) is the total opacity with ctt being 
the Thomson scattering cross-section and rup being the 
proton mass. 

For the equation of state we use pgas = (7 — l)e, where 
7 is the specific heat ratio. The temperature of the gas, 
T, can then be calculated from pgas = pk^T/prup, where 
fce is the Boltzmann constant and p is the mean molec¬ 
ular weight. To complete the set of equations, we apply 
flux limited diffusion (FLD) approximation developed by 
Levermore and Pomraning (1981). Adopting this ap¬ 
proximation, the radiation flux and the radiation pres¬ 
sure tensor are expressed in terms of the radiation energy 
density (Turner & Stone 2001, Ohsuga et al. 2005). 

In this paper, the dynamical viscosity coefficient is 
given by 

Pgas + AFq ^ ^ 

“ Ok ■ 

where a is the viscosity parameter, Ok is the Keplerian 
angular speed, and A is the flux limiter. In this form, the 
r(/?-component of the viscous stress tensor is 


(X 


aptotai (optically thick limit) 
apgas (optically thin limit) 


( 8 ) 


where ptotai is the total pressure. This is because the 
flux limiter. A, becomes 1/3 in the optically thick limit, 
and vanishes in the optically thin limit in the framework 
of the FLD approximation. The flux limiter. A, is al¬ 
most 1/3 in the disk region. Thus, our viscosity model is 
basically the same as the a prescription of the viscosity 
proposed by Shakura & Sunyaev (1973), although A is 
almost null above and below the disk. 



2D RHD MODEL FOR LIMIT-CYCLE DISK OSCILLATIONS 


3 


The computational domain consists of spherical shells 
of 3rg < r < 500rg and 0 < 0 < O.Stt, and is divided 
into 96 X 96 grid cells, where Vg is the Schwarzschild ra¬ 
dius. We start the calculations with a hot, rarefied, and 
optically-thin atmosphere. There is no cold dense disk 
initially, and we assume that matter is continuously in¬ 
jected into the computational domain through the outer- 
disk boundary (r = 500rg, 0 . 457 r < 0 < O.Stt). There¬ 
fore, we can avoid the influence of the initial configura¬ 
tion on the final result, although a long integration time 
is required. The injected matter is supposed to have a 
specific angular momentum corresponding to the Keple- 
rian angular momentum at r = lOOrg, and we set the 
injected mass-accretion rate (mass input rate) so as to 
be constant at the boundary. 

Throughout the present study, we assume the black- 
hole mass to be M = IOMq, a = 0.1, 7 = 5/3, and 
/r = 0.5. For the absorption opacity, we consider the free- 
free absorption and the bound-free absorption for solar 
metallicity (Hayashi, Hoshi, & Sugimito 1962, Rybicki & 
Lightman 1979). 

3. RESULTS 

Figure 1 represents the time evolution of the mass ac¬ 
cretion rate onto the black hole, Macc (blue), the out¬ 
flow rate, Mout (magenta), the luminosity, L (red), and 
the trapped luminosity, Ltrap (green). Here, the outflow 
rate means the ejected mass per unit time at the outer 
boundary of the computational domain. The luminos¬ 
ity is evaluated by integrating the radiation flux at the 
outer boundary. On the other hand, the trapped lumi¬ 
nosity is the integration of the radiation flux at the inner 
boundary, which means the radiation energy swallowed 
by the black hole per unit time. The luminosity, L, al¬ 
most equals to the luminosity of the disk, since the emis¬ 
sion of the less dense atmosphere surrounding the disk 
does not contribute very much to the luminosity. The 
mass input rate is set to be Minput = 100 Le/c^- 

The injected matter accumulates within the compu- 
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Fig. 1.— The time evolution of the mass accretion rate (blue), 
outflow rate (magenta), the luminosity (red), and trapped luminos¬ 
ity (green), for Minput = IOOLe/c^- Here, we employ the viscosity 
model, which gives the simple a prescription proposed by Shakura 
Sz Sunyaev (1973), tnp oc ptotah inside the disk. 


tational domain and the accretion disk forms around 
the black hole {t < 100 s). Subsequently, the limit- 
cycle oscillation starts. At f > 100s, the mass accre¬ 
tion rate onto the black hole (M^acc) drastically varies, 
although the mass is continuously injected into the 
computational domain at constant rate. It suddenly 
rises, retains high value for about 40 s, and then de¬ 
cays. Such time variation of the mass accretion rate 
causes the quasi-periodic luminosity changes between 
high- and low-luminosity states. Whereas the luminos¬ 
ity (L) is not more than 0.3Le in the low-luminosity 
state, it reaches around 2 Le in the high-luminosity 
state. The burst duration is also about 40 s, which 
is roughly correspond to the viscous timescale, tvis 
35 . 8 (M/lOM 0 )(r/lOOrg) 3 / 2 (a/O.l)-i(iL/O. 2 r )-2 s. The 
physical mechanism of such limit-cycle oscillations is the 
thermal instability in the radiation-pressure dominant re¬ 
gion (discussed later). 

Moreover, we can see in this figure that our simula¬ 
tions reveal two important phenomena. One of them 
is that the outburst is accompanied by the strong out¬ 
flow. As shown in the top panel, while the outflow rate 
(Afout) is very small in the low-luminosity state, it sud¬ 
denly rises together with mass accretion rate and retains 
the high value during the outburst. Since the luminosity 
is larger than the Eddington luminosity, the outflow is 
driven by the strong radiation force. Another one is that 
the photon-trapping effects appear during the burst. It 
is found that the trapped luminosity (Ltrap) is compara¬ 
ble to the luminosity in the high-luminosity state. Thus, 
considerable number of photons falls onto the black hole, 
reducing the apparent luminosity, although the disk is 
very luminous in this state. In contrast, the photon trap¬ 
ping does not work effectively in the low-luminosity state. 
The multi-dimensional effects are significant in the high- 
luminosity state. 

To understand the difference of the disk structure be¬ 
tween two states, we display the cross-sectional view of 
the density distribution in the low-(left panel) and high- 
luminosity states (right panel) in Figure 2. The velocity 
vectors are overlaid in the right panel. It is clear that 
the disk becomes geometrically thick and is accompanied 
by the outflow in the high-luminosity state, whereas the 
thin disk forms in the low-luminosity state. Further, the 
patchy structure as well as the circular motion appears 
within the disk in the high-luminosity state. They would 
be caused by the Kelvin-Helmholtz instability around the 
disk surface and by probably convection (see Ohsuga et 
al. 2005). It is stressed again that the multi-dimensional 
effects are significant in the high-luminosity state. 

In Figure 3, we plot the radiation temperature pro¬ 
files at the equatorial plane (top panel) and the effec¬ 
tive temperature profiles (bottom panel) in the high- and 
low-luminosity states. Here, the radiation temperature 
is evaluated as Tr = (Ao/a)^^"^, where a is the radiation 
constant, and it is roughly equivalent to the gas temper¬ 
ature around the equatorial plane. As shown in the top 
panel, the temperature of the disk is a few times larger 
in the high-luminosity state than in the low-luminosity 
state. We calculate the effective temperature by solv¬ 
ing the radiation transfer equation in the face-on view. 
We find that the effective temperature is roughly pro¬ 
portional to in the low state and in the high 
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Fig. 2.— The 2D density distributions in the low-(left) and high-luminosity states (right). The elapsed times are 166 s and 215 s, 
respectively. The arrows in the right panel indicate the velocity vectors. The adopted mass input rate and viscosity model are same as 
those in Figure 1. 


state except for the very vicinity of the black hole (see 
bottom panel). They are consistent with the relations 
for the standard and slim disks. 

Here, we need to remark that the sum of the mass ac¬ 
cretion rate and outflow rate is about half of the mass in¬ 
put rate on average. It means that the total mass within 
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Fig. 3.— The radiation temperature profiles at the equatorial 
plane (top) and the effective temperature profiles (bottom) in the 
high- and low-luminosity states. The adopted mass input rate and 
the viscosity model are the same as those in Figure 1. 


the computational domain monotonically increases with 
time. However, the increase of the mass mainly occurs 
in the outer part of the disk (r » lOOr^), and it does 
not influence the time evolution of the inner part of the 
disk (r < lOOrg) so much. Thus, we conclude that our 
simulations can reproduce the quasi-periodic oscillations 
of an unstable disk around the black hole. Much longer 
computation time is required, compared with the present 
simulations, for reproducing Mjnput = Vacc + Vout on 
average. 

Finally, we prove that the limit-cycle behavior in our 
simulations is triggered by the thermal instability in the 
radiation-pressure dominant region. According to the 
accretion disk theory, the limit-cycle oscillations occur 
when two conditions are satisfied. One of them is that the 
mass accretion rate is moderately high, Macc>TE/c^. 
Another one is that the viscous stress tensor is propor¬ 
tional to the total pressure. To investigate the physi¬ 
cal mechanism of the oscillations in our simulations, we 
study whether or not the limit-cycle oscillations occur 
by assuming Minput = 10 ^Le/c^ (comparison model A) 
and tr^p oc Pgas (comparison model B). The top panel 
in Figure 4 shows the time evolution of the mass ac¬ 
cretion rate as well as the luminosity of the comparison 
model A. The viscosity prescription is the same as that 
of the original model, that is, t^p oc ptotai in the disk. 
As shown in the panel, it is found that the disk becomes 
a quasi-steady (Mace 100Le/c^ and L ^ 3Le)- This 
represents that the disk is stabilized when the mass accre¬ 
tion rate highly exceeds the critical value. On the other 
hand, the disk with the viscosity model of oc Pgas 
does also not exhibit the limit-cycle oscillations. We plot 
the evolution of Mace and L of the comparison model 
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Fig. 4. — The time evolution of the mass accretion rate (blue) 
and the luminosity (red). In the top panel, the viscosity model is 
the same as that in Figure 1. But for high mass input rate, we 
adopt Minput = IO^Le/c^ (comparison model A). On the other 
hand, the mass input rate is the same as that in Figure 1. But for 
different viscosity law, we adopt tnp cx Pgas (comparison model B). 

B in the bottom panel. In this model, although the 
disk is not quasi-steady yet (Minput > Mace + M^out), 
the mass accretion rate is high enough to give rise to 
the instability in the case of tnp oc ptotai (Chen and 
Taam 1993, see also Kato et al. 1998), and the elapse 
time exceeds the viscous timescale in the inner region, 
Uis ~ 16O(M/lOM0)(r/5rg)3/2(a/O.l)-i(M/O.Olr)-2 s. 
It means that the viscosity model of t^p oc Pgas does not 
cause the disk instability. Since we succeed in reproduc¬ 
ing the prediction of the disk theory, we can conclude 
that recurrent outbursts in our simulations are caused 
by the disk instability in the radiation-pressure domi¬ 
nant region. 

4. DISCUSSION 

Resulting light curve in our simulations gives a nice fit 
to the time variation of the luminosity of GRS 1915-1-105. 
Based on the analysis of the data of GRS 1915-1-105 
taken by RXTE and ASCA, Yamaoka et al. (2001) have 
produced the light curve (see also Watarai & Mineshige 
2003). The luminosity is several times higher in the high- 
luminosity state than in the low-luminosity state. The 
duration of the high-luminosity state is about 30 — 50 s, 
and there is sharp edges between the high and low states. 
Our numerical results succeed in reproducing these ob¬ 
served features as shown in Figure 1. 

In our simulations, the observed bursting behavior is 
reproduced using the simple a prescription of the vis¬ 
cosity proposed by Shakura & Sunyaev (1973). [Note 
that our viscosity model gives trip oc ptotai inside the 
disk as we have already mentioned.] This also might be 
purely 2D effect, since the modification of the viscos¬ 
ity model is required by the ID simulations. Watarai 
& Mineshige (2003) explained the X-ray observations of 


GRS 1915-1-105 using Up oc ptotai(Pgas/ptotai)°'^- Some 
authors have proposed the phenomenological viscosity 
model, in which the effects of the dissipation of the en¬ 
ergy in disk corona or in jets are simply taken into con¬ 
sideration (Nayakshin et al. 2000; Janiuk et al. 2000, 
2002, 2005), to match the observed burst duration and 
luminosity. 

The multi-dimensional motion, which is treated in 
our simulations, would affect the time evolution of the 
disk. The mass, momentum/angular momentum, and 
energy are extracted from the disk by the outflow and 
transported by the convective motion and circulation. 
Through these multi-dimensional effects, the viscosity 
law might practically change and affect the limit-cycle 
behavior of the disk. To investigate the role of the multi¬ 
dimensional motion for the time evolution of the disk in 
detail is left as a future work. 

The radial profiles of the temperature and pressure 
of the disk in our (grid-based) simulations look similar 
to those of the SPH simulations (case 2 in Teresi et al. 
2004b). The temperature of the disk is a few times larger 
in the high-luminosity state than in the low-luminosity 
state (see top panel in Figure 3), and the radiation pres¬ 
sure exceeds the gas pressure by more than two orders 
of magnitude in the high state. However, the duration 
time of the low-luminosity state in their SPH simulations 
is quite differ from our results. In their simulations, the 
disk stays in the low state for a few second and it does not 
agree with observations of GRS 1915-1-105 (e.g. Janiuk 
& Czerny 2005). Our simulations indicate that the dura¬ 
tion time of the low state is roughly comparable to that 
of the high state. Also, the luminosity is a few times 
larger in their simulations than in our study, although 
the relatively small mass-input rate, ~ 30 Le/c^, is em¬ 
ployed in their simulations (we set Minput = 100 Le/c^ in 
the present work). These discrepancies may arise from 
the difference of the method of the numerical calcula¬ 
tions. Also, the treatment of the multi-dimensional ef¬ 
fects might lead to the difference of the time evolution 
of the disk. As we have already mentioned, the out¬ 
flow motion as well as the photon trapping are carefully 
treated in our simulations, although the SPH simulations 
by Teresi et al (2004b) focused on the time evolution of 
the disk only. 

Although we investigate the disk instability in the 
framework of the a prescription of the viscosity, the re¬ 
alistic formalism about the viscosity should be investi¬ 
gated from the magneto-hydrodynamical point of view 
(e.g., Machida et al. 2000; Stone & Pringle 2001). We 
need radiation MHD simulations as future work. 

5. CONCLUSIONS 

By performing the 2D RHD simulations, we investigate 
the time evolution of the accretion disks around the black 
holes and find the limit-cycle oscillations. In this study, 
we carefully treat the 2D effects including the outflow 
and the photon trapping. We summarize our results as 
follows: 

(1) The luminosity sharply rises from ~ 0.3Le to ^ 
2 Le in our simulations. Duration of the high-luminosity 
state is about 30 — 50 s. The resulting variation ampli¬ 
tude and duration nicely fit to the observations of micro¬ 
quasar, GRS 1915-1-105. 

(2) It is also found that the 2D effects are significant 
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in the high-luminosity state. In this state, the trapped 
luminosity is comparable to the luminosity due to the 
effective photon trapping. The outflow, circular motion, 
and patchy structure appear in this state. The outflow 
is driven by the strong radiation force. 

(3) The physical mechanism, which causes the limit- 
cycle oscillations, is the thermal instability in the 
radiation-pressure dominant region. In our simulations, 
the disk is stabilized when the mass accretion rate highly 
exceeds the critical value. The disk does not exhibit the 
oscillations if the viscous stress tensor is proportional to 


the gas pressure only. 
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